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ABSTRACT 


The  Single  Pass  Fit  Program  (SPFP)  is  designed  for  use  on  the  CDC  6400 
or  6600  machine  as  the  principal  computational  tool  for  preliminary  orbit 
determination  and  processing  of  single  arcs  of  raw  tracking  data.  By  an 
iterative  differential  correction  procedure,  it  solves  for  a  set  of  trajectory 
parameters,  namely,  initial  position  and  velocity  (and  drag  coefficient  if 
desired)  of  the  vehicle,  which  minimize  the  differences  between  measured 
and  computed  observation  measurements.  A  comprehensive  description  of 
the  mathematical  model  implemented  in  the  SPFP  is  presented  with  emphasis 
placed  on  key  characteristics,  such  as  reference  systems  and  equations, 
data  processing  and  self-initializing  techniques,  residual  editing,  and  differ¬ 
ential  correction.  In  addition,  a  basic  usage  guide  is  given  with  complete 
instructions  for  input  data  preparation  and  program  operation. 
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NOMENCLATURE 
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perturbative  acceleration  vector 

matrix  of  constraints 

set  of  refined  observations 

range  observation  measurement 

range  rate  observation  measurement 

unit  observation  vector 

magnitude  of  the  vehicle  velocity  vector 

station  frame  vector 

weighting  matrix 

principal  cartesian  (ECI)  axes 

semimajor  axis 

semiminor  axis 

eccentricity 

analytic  two- body  scalar  functions 

height 

inclination 


-IX- 


NOMENCLATURE  (Cont'd) 


N 

n 

o 

£ 

r 

£ 

t 

V 

x,y,z 

a 

P 

2 

6 

€ 

X 

P 

P 

cr 

T 


0' 


(y 


n 


total  number  of  observations 

number  of  accepted  observations 

observation  measurement 

vehicle  position  vector 

magnitude  of  the  vehicle  position  vector 

vehicle  velocity  vector 

time 

true  anomaly 

cartesian  (ECI)  coordinates 
right  ascension 
flight  path  angle 
vector  of  parameters 
declination 

small  positive  quantity 
geodetic  longitude 
gravitational  constant 
atmospheric  density 
weighting  factor 
time  of  perigee  passage 
geodetic  latitude 

argument  of  perigee,  or  rotation  rate 
right  ascension  of  the  ascending  node 


-x- 


1.  INTRODUCTION 


A.  PURPOSE,  SCOPE,  AND  LIMITATIONS 

This  document  is  intended  to  serve  as  .a  technical  reference  manual  and 
basic  usage  guide  for  the  Single  Pass  Fit  Program  (SPFP),  Preliminary 
orbit  determination,  as  defined  in  Section  H,  is  a  process  for  estimating 
from  tracking  data  the  initial -conditions  (position  and  velocity  of  a  vehicle) 
for  the  differential  equations  of  motion.  The  determination  process  imple¬ 
mented  in  the  SPFP  utilizes  a  precision  numerical  integration  method  to 
generate  the  vehicle's  trajectory  and  associated  partial  derivatives,  and  an 
estimation  technique  that  determines  the  initial  conditions  by  a  bounded  leasts 
squares  differential  correction  procedure.  A  comprehensive  description  of 
the  mathematical  model  used  is  presented  with  emphasis  placed  on  key  pro¬ 
gram  characteristics.. 

Sections  III  (Reference  Systems  and  Equations),  IV  (Functional  Compo¬ 
nents),  and  V  (Program  Usage)  are  sufficiently  detailed  to  meet  the  needs  of 
Aerospace  Corporation  technical  staff  members  concerned  with  preliminary 
orbit  determination..  In  the  interest  of  technical  completeness,  the  user  is 
referred  to  the  TRACE  Orbit  Determination  Program,  Version  D  (Rui.  1), 
which  covers  those  theoretical  concepts  omitted  in  this  document. 

B,  HISTORICAL  BACKGROUND 

In  the  past,  orbit  determination  has  been  performed  by  the  generalized 
TRACE -D  Orbit  Determination  Program  at  Aerospace  Corporation.  Tice 
SPFP  was  developed  because  of  the  increasing  need  for  a  specialized,  com¬ 
pact,  and  rapid  computational  tool  to  accommodate  preliminary  orbit  deter¬ 
mination  and  processing  of  single  arcs  of  raw  tracking  data.  Its  design 
char  cteristlcs  (written  entirely  in  6600  FORTRAN)  closely  parallel  those  of 
TRACE-66,  the  latest  version  of  TRACE. 


c. 


SUMMARY  OF  KEY  FUNCTIONS 


By  an  iterative  differential  correction  procedure,  the  SPFP  solves  for 
a  set  of  parameters,  namely,  position  and  velocity  (and  drag  coefficient  if 
desired)  of  the  vehicle  that  minimize  the  differences  between  measured  and 
computed  tracking  observations  from  single  arcs  of  data.  The  preliminary 
orbit  determination  process  used  in  SPFP  can  be  identified  by  the  following 
functions : 

1.  Data  processing  refines  the  raw  observation  measurements  by 
conversion  and  pre-editing  (or  rejection)  criteria.  Range  and 
elevation  measurements  are  subject  to  maximum  and  minimum 
constraints,  and  scaling  and  refraction  corrections. 

2.  Self-initializing  of  the  initial  position  and  velocity  from  selected 
data  points  is  performed  with  a  gaussian  iterative  scheme. 

3.  Trajectory  generation  of  position  and  velocity,  and  associated 
partial  derivatives  at  each  data  point,  is  accomplished  by  numer¬ 
ical  integration  of  the  equations  of  motion  in  a  Cowell  formulation 
with  an  eighth-order  differencing  method. 

4.  Residual  evaluation  (that  is,  measurement  residuals)  is  made 
from  the  differences  between  observed  and  computed  values. 

5.  Observation  partial  derivative  evaluation  is  the  computation  of 
the  derivative  of  each  data  measurement  with  respect  to  trajec¬ 
tory  parameters. 

&<  Residual  editing  is  performed  on  those  points  v/hose  residuals 

exceed  a  predetermined  multiple  of  the  prior  iterative  rms  (roc.t 
mean  square)  value. 

7.  Matrix  accumulation  consists  of  storing  observation  partials  and 
residuals  in  a  matrix  A,  which  is  used  to  form  the  weighted 
normalized  symmetric  matrix  product  A^A. 

8.  Differential  correction  of  parameters  depends  upon  a  bounded 
least-squares  process  that  minimizes  the  sum-of-squares  of  the 
observation  residuals. 


SECTION  II 


PRELIMINARY  ORBIT  DETERMINATION 

A.  THE  BASIC  PROBLEM  AND  KEY  CHARACTERISTICS 

OF  ORBiT  DETERMINATION 

Fundamentally,  orbit  determination  is  a  process  of  statistically 
estimating  from  tracking  data  the  initial  state  vector  of  a  vehicle,  that  is 
of  determining  the  initial  conditions  (position  and  velocity  vectors)  for  the 
vector  differential  equation  of  motion.  The  basic  problem  of  preliminary 
orbit  determination  is  to  establish  a  process  for  formulating  an  approxima¬ 
tion  of  the  initial  state  vector.  Many  such  processes  exist  at  various  levels 
of  complexity.  In  the  case  of  the  Single  Pass  Fit  Program,  the  preliminary 
determination  of  satellite  orbits  from  single  arcs  of  data  is  accomplished 
with  a  weighted  least- squares  differential  correction  procedure. 

Orbit  determination  is  characterized  by  the  trajectory  and  partial 
derivative  generation  methods,  and  by  the  mathematical  scheme  that 
estimates  the  initial  state  vector  of  the  vehicle  from  differences  between 
the  observed  and  computed  measurements.  For  preliminary  orbit  deter¬ 
mination,  the  SPFP  utilizes  a  precision  numerical  integration  method  to 
generate  the  trajectory  and  associated  partial  derivatives;  estimation  of 
the  initial  state  vector  is  performed  by  a  bounded  least- squares  differential 
correction  procedure  that  uses  residual  editing  techniques  to  enhance 
convergence. 

B.  THE  PRELIMINARY  DETERMINATION  MODEL 

The  computational  model  used  in  SPFP  to  determine  orbits  from 
tracking  data  can  be  characterized  by  four  specific  elements:  simulated 
vehicle  environment,  parameter  specification,  data  and  residual  editing, 
and  weights. 
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SIMULATED  VEHICLE  ENVIRONMENT 


To  simulate  the  environment  in  which  a  vehicle  is  moving,  the  following 
physical  forces  are  considered  to  be  acting: 

o  Gravitational  attraction  of  an  aspherical  earth  derived  from  a 
generalized  potential  function. 

e  Atmospheric  drag  (U.  S.  Standard  Atmosphere  1962,  COESA, 
density  model,  Ref.  2). 

Inclusion  of  any  subset  of  these  forces  is  completely  at  the  discretion 
of  the  program  user  through  input  options. 

2.  PARAMETER  SPECIFICATION 

The  SPFP  is  specifically  designed  to  determine  by  differential  correc¬ 
tion  the  initial  state  vector  Sq  (r,  r)  and,  if  desired,  the  drag  coefficient 

cda/w. 

3.  DATA  AND  RESIDUAL  EDITING 

The  following  are  raw  earth-fixed  data  observations  accepted  by  SPFP: 

9  Slant  range  (R) 

•  Local  azimuth  (A) 

e  Local  elevation  (E) 

e  Range  rate  (R) 

These  are  subject  to  preprocessing  editing  criteria,  such  as  maximum- 
\minimum  range  and  elevation,  and  data  arc  (pass)  duration.  Also,  each 
time  a  differential  correction  cycle  is  performed,  the  residuals  are  subject 
to  editing  criteria  determined  by  the  root- mean-square  (rms)  value  for  all 
residuals  in  the  previous  iteration. 

4.  WEIGHTS 

In  general,  measurements  are  assigned  a  priori  weights  according  to 
tracking- station  and  observation- measurement  type.  The  twofold  purpose  of 
the  weights  is  to  normalize  units  so  that  different  types  of  measurements  can 
be  mixed  in  the  least  squares  process,  and  to  adjust  their  relative  influence 
in  fitting  measurements  of  different  accuracies. 
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SECTION  in 


REFERENCE  SYSTEMS  AND  EQUATIONS 


The  basic  reference  systems  and  equations  used  in  the  computational 
method  of  the  SPFP  are  now  defined  and  discussed,  preceded  by  the  symbols 
used  throughout  this  report. 

A.  SYMBOLS 

.  time  differentiation  and  vector  dot  production  operation 

•  •  time  differentiation 

a  vector;  for  example,  £  is  the  vehicle  position  vector 

(  )  relates  to  a  functional  relationship,  such  as  S(r,  r)  the  vehicle 
state  vector,  which  is  a  function  of  £  and  £,  or  S(a,  e,  i,  £2,  co, 
t)  the  state  vector  expressed  in  terms  of  elliptical  elements. 

Also,  used  to  denote  row  vectors;  for  example,  £  =  (x,  y,  z) 

~  raw  observation  or  measurement 

A  increment  or  difference 

x  vector  cross  product  operation 

|  |  absolute  value,  or  vector  magnitude  operation;  for  example,  r  =  |£  [ 
V  vector  gradient  operator 

/\j  r\j  «\»  *\j 

{  }  a  matrix,  or  a  set  of  elements;  for  example,  Oj  =  {  Tj,  Rj,  Aj,  E^, 
the  set  of  raw  measurements  of  the  i^  observation 

j|  |J  the  vector  norm  operator 

l[  ]  the  integer  operator  of  a  scalar 
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a  matrix;  for  example. 


[  ] 


r  =  [x,  y.,  z]  (row  vector)  or  r  = 


x 

y 

z 


(column  vector) 


I  = 


1  0  0 

0  1  0 

.  0  0  1 . 


3  by  3  identify  matrix 


A  =  [  A]  4  by  7  matrix  of  partials  and  residuals 
Superscripts 

*  particular  value  of  a  predetermined  constant 

s  the  station 

i  particular  index  relative  to  a  set  or  vector 

c  computed  or  simulated  value 

T  matrix  transpose  symbol 

Subscripts 

o  initial  value  or  reference  point 

g  the  Greenwich  meridian 

s  the  station 

i,  j,k,  l  particular  index  denoting  members  of  a  set  or  components 
of  a  vector.  In  general,  i,  j,  k,  i  =  i,  2,  .  .  . 

e  the  earth 

a  a  transformation  that  is  a  function  of  a 

$  a  transformation  that  is  a  function  of  $ 

y  differentiation  of  a  vector  with  respect  to  the  parameter  y 
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B.  REFERENCE  SYSTEMS 


Preliminary  orbit  determination  is  a  computational  process  of  coordinate 
transformations  and  requires  selection  of  the  proper  reference  systems  for 
the  performance  of  the  necessary  function.  Several  such  systems  are  used 
in  the  SPFP;  basic  to  all  is  the  earth-centered  inertial  (ECI)  system,  which 
has  as  its  fundamental  plane  and  principal  axis  the  true  equator  at  epoch  and 
;he  mean  equinox  at  midnight  of  the  date  of  epoch.  This  reference  system 
does  not  account  for  precession  and  nutation  effects. 

1.  EARTH-CENTERED  INERTIAL  SYSTEM 

Inertial  position,  velocity,  and  acceleration  vectors  are  oriented  in 
the  reference  system  illustrated  in  Figure  i. 

The  position  and  velocity  components  of  a  vehicle's  state  vector  within 
this  reference  system  may  in  addition  be  expressed  in  terms  of  the  spherical 
or  classical  elliptic  systems. 


Figpre  1.  Earth-centered  Inertial  Reference  System 
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In  Figure  i: 

O  is  the  origin,  which  coincides  with  the  center  of  the  earth 
r_  is  the  radius  vector  (x,  y,  z) 
is  the  velocity  vector  (x,  y,  z) 

X  is  the  direction  of  the  vernal  equinox,  T  in  the  equator 
plane 

Y  is  the  right-hand  axis  to  X  and  Z 

Z  is  the  direction  north  perpendicular  to  the  equator  plane 
x  is  the  position  component  in  the  X  direction 

y  is  the  position  component  in  the  Y  direction 

z  is  the  position  component  in  the  Z  direction 

x  is  the  velocity  component  in  the  X  direction 

y  is  the  velocity  component  in  the  Y  direction 

z  is  the  velocity  component  in  the  Z  direction 

2.  SPHERICAL  SYSTEM 

In  the  spherical  system  shown  in  Figure  2,  the  vehicle  state  vector  at 
time  t  may  be  written  as 


S  =  S(a,  5,(3,  A,  r,  Y) 


Figure  2.  Spherical  Coordinate  Refercr.*.e  System 
In  Figure  2: 

a  is  the  right  ascension 

6  is  the  geocentric  latitude 

P  is  the  angle  between  radius  and  velocity  vectors 

A  is  the  azimuth  angle  between  the  projection  of  the  velocity 
vector  in  the  plane  normal  to  the  geocentric  radius  vector 
and  the  direction  to  due  north 

r  is  the  magnitude  of  the  geocentric  radius  vector  _r 

V  is  the  magnitude  of  the  velocity  vector  i_ 

a  is  the  right  ascension  of  Greenwich 

g  6 

X.  is  the  earth-fixed  longitude  of  the  vehicle 
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3. 


CLASSICAL  ELLIPTIC  SYSTEM 


The  state  vector  in  the  classical  elliptic  system  (Fig.  3)  is  given  at 
time  t  by 


S  =  S(a,  e,  i,f2,  u,  t) 


Figure  3.  Classical  Elliptic  Reference  Frame 
In  Figure  3: 

a  is  the  semimajor  axis  of  the  osculating  conic 
e  is  the  eccentricity- 

i  is  the  inclination 

Q  is  the  right  ascension  of  the  ascending  node 

w  is  the  angle  between  the  direction  of  perigee  and  the  line 

of  nodes 


r  is  the  time  of  last  perigee  passage  referenced  to  midnight 
of  epoch  date 

v  is  the  true  anomaly 

b  is  the  semiminor  axis 

4.  STATION  SYSTEM  (W -FRAME) 

The  station-dependent  measurements  are  transformed  to  this  system 
via  an  intermediate  reference  frame  (Fig.  4)  defined  by  the  station  coordin¬ 
ates  (X  ,  <p' ,  h).  If  Ws  denotes  the  station,  then  the  W -frame  is  the  plane 
s 

containing  the  vector  Ws  =  (W^,  so  that  =  0  and  coincides 

with  the  Z  axis. 


Figure  4.  Station  Reference  Frame 


In  Figure  4: 

\g  is  the  longitude  of  the  station 

<f>'  is  the  geodetic  latitude  of  the  station 

h  is  the  height  of  the  station  above  the  earth  ellipsoid 
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5. 


OBSERVATION  FRAME 


The  measurements  of  the  vehicle's  position  with  respect  to  the  station 
are  made  in  the  frame  defined  by  the  local  horizontal  plane  that  is  tangent  to 
the  earth  ellipsoid  (fig.  5). 


Figure  5.  The  Observation  Frame 


In  Figure  5: 

R  is  the  distance  between  the  station  and  the  vehicle  (range) 

A  is  the  angle  between  the  projection  of  the  range  in  the  local 

horizon  plane  and  due  north 

E  is  the  angle  between  the  projection  and  the  range  direction 
Ws  is  the  origin  of  the  observation  frame  and  is  related  to  the 

g 

station  frame  by  the  W  vector 

C.  TRANSFORMATION  EQUATIONS 

The  transformation  of  the  state  vector  to  the  other  systems  involves 
the  following  equations. 
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L  FROM  ECI  TO  THE  SPHERICAL  SYSTEM 


S(x,y,ZjX,y,z)  to  S(a4  6,  |3,  A,  r,  V) 

a  =  arc  tan  (y/x) 

5  =  arc  sin  (z/r) 

p  =  arc  cos  [  (xx  4  yy  4  zz)/rV] 

A  -  arc  tan  jr(xy  -  yx)/[y(yz  -  zy)  -  x(zx  -  xz)]J 

f  2  ,  2  2,2 

r  =  (x  4  y  4  z  } 

/  -  2  ^  2  .2,2  /ii 

Y  =  (x  4  y  4  z  )  (1) 

2.  FROM  ECI  TO  THE  CLASSICAL  ELLIPTIC  SYSTEM 
S(x,  y,  z,  x,  y,  z)  to  S(a,  e,  i,S2,  uij  -r) 

a  =  (2/r  -  Y2/^)"1 

1^ 

e  =  [  (e  cos  E)2  4  (e  sin  E)2]2 

A 

i  =  arc  tan  [(pj  +  -  P^l] 

S2  =  arc  tan  [(^2^3  “  ^*3^2 ^^1^3  ”  ^3®P" 
w  =  arc  tan  (P^/Qg) 

t  =  t  -  M/n  (2) 
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where 


y 

e  cos  E 


,2  2  2,2 

(x  +  y  +  z  ) 

.2.2.2 

x  +  y  +  z 


‘-f 


2 

e  sin  E  =  (xx  +  yy  +  zz)/(  |  a  |  p) 

p  =  [  rZY2  -  (xx  +  yy  +  zz)2]  /\i 

D  =  (xx  +  yy  +  zz)/ep. 

D  -  e  cos  E/er 

H  =  (r  -  pj/ex/ilp 

H  =  (xx  +  yy  +  zz)/re  s/pip 

P.  =  Dx  -  Dx 


P 

3 

Qi 

°2 

°3 

n 

M 

E 


=  Dy  -  Dy 
=  Dz  -  Dz 
=  Hx  -  Hx 
=  Hy  -  Hy 
=  Hz  -  Hz 

J. 

=  ((*/|a|3)2 

=  E  -  e  sin  E 

=  arc  tan  (e  sin  E/e  cos  E) 
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3.  THE  STATION  COORDINATES  TO  THE  W- FRAME 


(\g50',h)  to  (W^W^W^) 


W ®  =  (a-gae  +  h)ccs 

w2s0 

W3  =[(i  -  f)2bgae +  h]sin^f  (3) 


where 

ag  is  the  mean  equatorial  radius  of  the  earth 

f  is  the  earth  flattening  coefficient 

-2  2.  ..  -.2  .  2 
a  =  cos  <j>  +  (1  -  f)  sin  0 
s 

i-2  .  2  ,  ..  ..2  2  , 

b  =  sin  0  +  (1  -  f)  cos  0 

s 

4.  OBSERVATION  MEASUREMENTS  TO  THE  ECI  SYSTEM 
(R,  A,  E)  to  (x,  y,  z) 


The  transformation  of  observation  measurements  to  the  ECI  system 
permits  the  computation  of  the  initial  state  vector  by  a  gaussian  two-point 
scheme  with  an  analytic  two-body  formulation.  Thus 


(cos  a  -sin  a  0^ 
sin  a  cos  a  0 


\  0 


0  1 


r 

0 

-sin  0 ' 

COS  0'^ 

(r  cos  E 

sin  A^ 

K\ 

1 

0 

0 

R  cos  E 

cos  A 

W2 

.\° 

COS  0' 

sin  <p'  1 

^R  sin  E 

i 

K/ 

(4) 
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where 


a  =  a  +  to  (t  -  t  ) 

gQ  s  e'  o' 

where 

a  =  the  right  ascension  of  Greenwich  at  midnight  of  epoch  date 
8o 

w  =  the  earth  rotation  rate 
e 

t-tQ  =  the  elapsed  time  since  midnight  of  epoch  date 
5.  THE  STATE  VECTOR  TO  THE  W-FRAME 
(x,  y,  z)  to  (-W  W z>  W3) 

To  compute  simulated  measurements,  the  state  vector  is  transformed 
to  the  W -frame  by 
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D.  TRAJECTORY  EQUATIONS 

Three  sets  of  trajectory  equ^aons  are  used  in  the  SPFP  to  compute  the 
vehicle  state  vector,  S(r,  or  associated  partial  derivatives  of  the  state  at 
some  time  t.  T hc«e  equation  sets  are  described  by  their  functional  usage  in 
Table  I. 


Table  I.  Functional  Use  of  Equation  Sets 

Equation  Set 

Usage  Function 

Two -body  Equations 

To  evaluate  the  state  vector  SQ(r,  r) 
at  some  specified  epoch  time  t  from 
initial  conditions  given  at  observation 
time  t^ 

Equations  of  Motion 

To  generate  the  state  vector  S(r,jr) 

at  each  observation  time  t. 

l 

Variational  Equations 

To  generate  partial  derivatives  of  the 

state  vector  with  respect  to  the  param¬ 
eter  vector  at  each  observation  time 

t. 

l 

TWO-BODY  EQUATIONS 


Exact  closed-form  expressions  for  the  f  and  g  series  (Ref.  3)  are 
utilized  to  evaluate  the  initial  state  vector  SQ(r,  r)  at  some  specified  time  tQ, 
if  the  position  and  velocity  vectors  at  some  time  t^  are  given.  Let  r^  and 
r^  denote  the  given  vector  at  t^.;  then  the  initial  state  vector  at  t  can  be 
expressed  as 


io  =  +  s4  <7> 

and 

4  *  +  *4  (8) 
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where  f,  g,  f,  and  g  are  scalar  coefficients  given  by 

■) cos  ®  -  (tt)  V“ 


.  (-*' 


V 


sin  $ 


g 


■(f) 


r  sin  $ 
o 


$  =  (v  -  vk)  =  the  time  anomaly  difference 

(vq  -  v^)  =  function  of  [  (Eq  -  E^)]  =  the  eccentric  anomaly  difference 


n(to  -  VJ 

1  =  (E, 

-  Ek> +  \ 

o 

Jl. 

■Jfc  '  -k^ 

CV 

W 

V-  ) 

/u.  \ 

1/2 

n 

=  fe 

1 

'E  -  E, 
o  k 


sin(E  -  E.  ) 
'  o  k7 


M- 


=  product  of  the  Newtonian  gravitational  constant  and  mass 
of  the  earth 


g 


$ 


00k(*°si 


(r  cos  $  -  r  sin  <£  $) 
o  o 


sin  $  +  rQcos  $  <j?) 


sin  $  +  r  cos  $  $) 


h 

T 

r 

o 


l^k  X^kl 
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2. 


EQUATIONS  OF  MOTION 


The  state  vector  S(r,  r)  is  generated  at  each  observation  time  by  numeri¬ 
cally  integrating  the  equations  of  motion  in  a  Cowell  formulation  with  time  as 
the  independent  variable.  The  equations  of  motion  are  expressed  by  the 
second-order  vector  differential  equation 


(9) 


with  the  initial  values 


r(t  )  =  r 
— '  o  — o 


f(t  )  =  f 
—  o  — o 


where 

£  is  a  3-vector  of  ECI  components  of  position  (x,  y,  z) 

r  is  the  magnitude  of  r_  (r  =  |  r  | ) 

|j.  is  the  product  (GM)  of  the  Newtonian  gravitational  constant  and 
mass  of  the  earth 

F  is  the  perturbative  acceleration 

3 

In  particular,  the  primary  term  -  pr/r  is  the  inverse  square  central 
force  due  to  gravitation.  The  perturbative  term  F  is  the  sum  of  the  gravi¬ 
tational  acceleration  due  to  the  noncentral  force  field  of  the  earth  F^  and  the 
acceleration  due  to  atmospheric  drag  Fg. 

F.  is  derived  from  the  generalized  potential  function 


U  *  f  [‘  +  X)(t)  X)  P”  (sin  *)(Cnmcos  mX  +  Snm8in  mX)]  (10) 
n=2  m=0 


-19- 


where 


r. 


e 

3m 

n 


is  the  product  of  the  Newtonian  gravitational  constant 
and  the  earth  mass 


r,  <j>,  X  are  the  geocentric  distance,  latitude,  and  (east) 
longitude  of  the  vehicle 

is  the  mean  equatorial  radius  of  the  earth 


is  the  Legendre  associated  function  of  the  first  kind 
of  degree  n  and  order  m 


C  ,S 
nm  nm 


are  numerical  coefficients 
The  vector  form  of  FA  can  be  expressed  as 


Fa  =  VU(r)  +^= 
r 


(11) 


andFg,  the  force  due  to  atmospheric  drag,  is  expressed  by  the  vector 
equation 


Ib  =  -p 


fr)  (¥) 


(12) 


where 


p  is  the  density  at  height  h  above  an  oblate  earth 

CjjA/W  is  the  drag  coefficient 


is  the  vehicle  velocity  vector  relative  to  a  rotating 
atmosphere  [i_A  =  C*A,yA,  z A] ) 

=  £a  1 
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where 


iA  =  *  +  <v 

yA  =  y  -  “ax 

•  • 

z  .  =  z 

A 

u)^  =  rotation  rate  of  the  atmosphere 

The  density  model  used  to  obtain  p  in  the  SPFP  is  the  1962  U.  S. 
Standard  Atmosphere  (Ref.  2),  which  is  designed  to  represent  atmospheric 
density  under  average  conditions  of  solar  activity  through  middle-latitude 
bands.  It  is  similar  to  the  ARDC  1959  model  (Ref.  4)  in  which  the  sequence 
of  connected  linear  segments  involving  variation  of  molecular  scale  temper¬ 
ature  with  altitude  is  considered.  The  1962  model  reflects  density  values 
somewhat  lower  than  the  ARDC  1959  at  altitudes  above  200  km. 

3.  VARIATIONAL  EQUATIONS 

To  generate  partial  derivatives  oi  the  state  vector  with  respect  to  a 
parameter  set  jy,  variational  equations  are  numerically  integrated  and  then 
the  partial  derivatives  are  interpolated  at  each  observation  time.  In  the 
SPFP  only  parameters  associated  with  the  equations  of  motion  are  consid¬ 
ered;  in  particular,  initial  condition  parameters,  such  as  SQ(r, _r)  and 
Cj^A/W,  are  the  only  ones  specified  in  the  differential  correction  process. 

If  the  parameter  set  is  denoted  by  the  vector  it  can  be  functionally 
related  to  the  vector  equation  of  motion  by  the  form  of  Eq.  (9)  as 

r  =  r(r,r,^,  t) 
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Let  rY 


denote  9r/3y;  then  the  general  vector  variational  equation  is 


"Ml  *  j£ 
3rJ  — ' y  +  3y 


(13) 


with  the  necessary  initial  conditions 


i/V  *  <14> 

depending  on  the  parameter  y  of  X'  The  9F/3y  in  Eq.  (13)  is  called  the  non- 
homogeneous  term  of  the  variational  equation  for  y  of  Note  that  the  con¬ 
tents  of  the  brackets  in  Eq.  13  are  3  by  3  matrices.  The  matrix  in  the  first 
set  of  brackets  is  evaluated  as  the  sum  of  [v]  +  [T],  where 


(15) 


and 


(T]  = 


3r 


(16) 


represent  the  dependence  of  the  gravitational  and  atmospheric  drag  accelerations 
upon  the  position  of  the  vehicle.  The  matrix  [  3F/3r]  can  be  written  as 


r  r 
-A-A 


L  A 


(17) 


Initial  Conditions  for  Variational  Equations 


Initial  conditions  for  the  second- order  linear  vector  differential  given 
in  Eq.  (14)  are 


m,  •[£] 


=  I  (3  by  3  identify  matrix) 


\W~\  =  ^  (zero  matrix) 

L  — h oJt  L  — oJt 


where  %  =  (r^j^),  representing  SQ(r,r).  For  y  of  equal  to  (C^A/W) 


3CdA7W  t 


dL 

=  9C  A/W  .  = 

o  u  o 

L  j 


b.  Nonhomogeneous  Term  for  Variational  Equation 


The  nonhomogeneous  term  9F/9y  *n  Eq.  (13)  is  equal  to  zero,  if 
Y  =  xq,  yQ,  zQ,  xq,  y  ,  or  zq;  and  equals  F^CpA/W)"1  only  if 

y  -  (CjjA/W),  a  non-zero  quantity. 

E.  OBSERVATION  MEASUREMENTS  AND  ASSOCIATED  PARTIAL 


DERIVATIVES 

Equations  for  the  computation  of  simulated  measurements  and  their 

s 

associated  partial  derivatives  are  formulated  directly  from  the  W  vector. 
SIMULATED  OR  COMPUTED  MEASUREMENTS 


Range  (R)  =  *  Q3) 

Azimuth  (A)  =  arc  tan.(vi /v?) 


[•)« 


Elevation  (E) 


ASSOCIATED  PARTIAL  DERIVATIVES  OF  THE  MEASUREMENT  WITH 
RESPECT  TO  THE  PARAMETER  VECTOR 


Range 


OR 

9V 


Azimuth 


9A 

9y 


(confc) 


SECTION  IV 


FUNCTIONAL  COMPONENTS 


Seven  basic  functional  components  identify  the  computational  model  used 
in  the  SPFP  to  perform  preliminary  orbit  determination.  Figure  6  schemat¬ 
ically  illustrates  the  function  of  each  component  in  relation  to  the  solution  of 
the  preliminary  orbit  reconstruction  problem.  These  components  are  defined 
by  the  key  functions  they  perform,  as  shown  in  Table  II. 

Table  II.  Key  Functions  of  Basic  SPFP  Components 


COMPONENT  KEY  FUNCTION 


Input  Data  Processor 

Process  input  data 

Data  Processor 

Convert  and  refine  raw 
observation  measurements 

Self -Initializer  for  the  State  Vector 

Estimate  the  initial  state 
vector  from  selected  data 
points 

Trajectory  Generator 

Produce  the  state  vector  and 
its  associated  partial  derivatives 
at  each  data  point 

Residual  and  Partial  Derivative  Evaluator 

Evaluate  data  residuals  and 
data  partial  derivatives 

Residual  Editor  and  Matrix  Accumulator 

Edit  data  based  on  residual 
level,  and  accumulate  matrix 
of  normalized  ^partial 
derivatives  A^A 

Differential  Corrector 

Perform  least  square  calcula¬ 
tions  to  compute  the  parameter 
correction  vector 

TRAJECTORY 

GENERATOR 

LOOP 


RESIDUAL  AND 
PARTIAL  DERIVATIVE 
EVALUATOR 


ITERATIVE 
LOOP 
OF  THE 
DIFFERENTIAL 
CORRECTION 
PROCESS 


RESIDUAL  EDITOR 
AND  MATRIX 
ACCUMULATOR 


DIFFERENTIAL 

CORRECTOR 
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A. 


INPUT  DATA  PROCESSOR 


The  basic  function  of  the  input  data  processor  is  to  initialize  the 
observation-dependent  quantities  and  model-dependent  variables. 

i.  OBSERVATION-DEPENDENT  QUANTITIES 

The  observation-dependent  quantities  of  the  SPFP  consist  of  station 
location  and  observation  information,  editing  criteria,  iteration  and  initial 
condition  selection  controls,  and  conversion  factors. 

Station  location  and  observation  information  is  initialized  for  each  single 
arc  of  tracking  data.  This  'consists  of  items  such  as: 

©  Earth-fixed  station  location  coordinates 

©  Measurement  reference  data. 

©  Observation  data  format  and  measurement  cedes. 

Editing  criteria  include  quantities  such  as 

e  Range  and  elevation  refraction  corrections 

•  Range  scaling  factor 

•  Maximum -minimum  range  and  elevation  constraints 
©  Initial  residual  limits  for  each  measurement. 

Iteration  and  initial  condition  selection  controls  consist  of: 

«  Maximum  number  of  iterations 

•  Data  point  selection  index  for  self-determination 
of  initial  state  vector 

©  Measurement  weighting  sigmas 

©  Bounds  for  corrections  to  the  state  vector  estimate 

•  Scaling  factors  that  reduce  (increase)  bounds  on  a  diverging 
(converging)  iteration 

©  Relative  convergence  criterion. 

Conversion  factors,  such  as  distance  in  feet  per  earth  radii  (ft/er) 
and  velocity  (ft/sec/er /min)  constants,  are  necessary  for  input-output. 
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2. 


MODEL -DEPENDENT  VARIABLES 


The  model -dependent  variables  of  the  SPFP  consist  of  physical  constants- 
numerical  integration  controls,  and  force  model  specifications. 

Physical  constants  —  geopotential  and  atmospheric  model  constants  —  that 
define  the  physical  environment  of  the  vehicle  are: 

«  Gravitational  constant  (p) 

t>  Numerical  seopotential  coefficients  (C  ,  S  ) 

•  Mean  equatorial  radius  (ae) 

e  Earth  flattening  term  (f) 

o  Atmospheric  'constants . 

Numerical  integration  controls  consist  of: 

e  Initial  integration  step  size 

©  Maximum  and  minimum  allowable  step  size 

•  Ratio  of  Runge-Rutta-to -Cowell  step  size 

The  force  modal  is  further  specified  by  the  user  through  data  inputs  such  as: 

o  Number  of  terms  in  geopotential  model 
©  Geopotential  normalization 

o  Drag  coefficient  (CjyA/W) 

B.  DATA  PROCESSOR 

Raw  observation  measurements  are  subject  to  conversion  and  pre¬ 
editing  (or  rejection)  criteria  through  the  SPFP  data  processor.  The 
procedure  used  refines  the  raw  observation  set  Ointo  an  observation  set  O, 
which  is  directly  utilized  in  the  reconstruction  scheme. 

1.  RAW  OBSERVATION  SET,  6 


where 


N  is  the  total  number  of  observations  for  a  given  arc  of  data 

til 

T.  is  the  date  and  time  of  day  of  the  i  observation;  that  is, 

1  T.  =  {YEAR,  MONTH,  DAY,  HOUR,  MINUTE,  SECOND}. 

R^  is  the  observed  range  (ft) 

A^  is  the  observed  azimuth  (deg) 

E.  is  the  observed  elevation  (deg) 

R^  is  the  observed  range  rate  (ft/sec) 

til 

i  is  the  subscript  indicating  the  i  observation. 

2.  CONVERSION 
th 

Time  of  the  i  raw  observation  T^  is  the  first  quantity  to  be  refined. 
A  reference  Julian  date  (J.  D. )  is  computed  from  T.  or  is  prespecified  in 
order  to  redefine  each  Tj  into  t.  with  units  of  minutes  from  midnight  of  the 
reference  Julian  date. 

Observed  measurements  have  the  conversion  units  shown  in  Table  HI. 


Table  HI.  Observed  Measurement  Conversion  Units 


MEASUREMENT 

RAW  UNIT 

REFINED  UNIT 

Range 

feet 

earth  radii 

Azimuth 

degrees 

radians 

Elevation 

degrees 

radians 

Range  Rate 

feet/sec 

earth  radii /minute 

3.  PRE-EDITING 


In  the  usual  reconstruction  situation,  range  and  elevation  measurements 
are  required  to  lie  within  specified  limits.  Thus  each  R.  and  E.  are  sub-  ; 
jected  to  the  following  predetermined  constraints  for  i  =  1,  .  .  .  ,  N; 
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rmin 

rmax 

emin 


Minimum  range  constraint 
Maximum  range  constraint 
Minimum  elevation  constraint 


That  is,  any  raw  range  or  elevation  measurement  not  satisfying  the  above 

constraints  deletes  its  corresponding  observation  from  6;  hence,  N  is 

reduced  by  one  for  each  rejected  observation. 

The  data  arc  is  sorted  and  examined  to  determine  if  At.  =  t.  -  t, 

1  1  1  * 

(i  =  1,  .  .  .  ,  N)  is  greater  than  some  prespecified  arc  duration  value.  At  . 

Sic 

If  At.  >  At  ,  then  the  i  observation  is  deleted  and  N  is  reduced  by  one. 

4.  REFINED  OBSERVATION  SET  O 


After  the  measurements  for  each  observation  are  converted  and  pre¬ 
edited,  the  range  data  is  subject  to  scaling  and  refraction  corrections  upon 
option;  in  addition,  a  refraction  correction  may  be  applied  to  the  elevation 
measurement.  Therefore,  upon  option 

E.  =  E.  +  AE  ,  if  E.  S  C.  1  radian 

li  i 

E.  =  E.  +  AE^  ,  if  E.  >  0. 1  radian 

ii  l 

jj; 

where  AE  and  AE  denote  predetermined  refraction  corrections. 

The  resultant  set  of  observations  is  defined  as  the  refined  observation 
set  O;  that  is 


O  = 


i  —  1,  ...,  n 


(21) 


where 

n  is  the  total  number  of  acceptable  observations 

th 

t.  is  the  time  of  the  i  data  point  in  minutes 

1  from  midnight  of  the  reference  Julian  date 

R.  is  the  range  (er) 
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1 


A.  is  the  azimuth  (rad) 

E^  is  the  elevation  (rad) 

R.  is  the  range  rate  (er/min) 

th 

i  is  the  subscript  indicating  the  i  observation 

with  the  set  Ch  =  {t^,  R^,  A.,  E^,  R.},  denoting  the  i  data  point. 

C.  SELF-INITIALIZATION  OF  THE  STATE  VECTOR 

The  technique  used  in  the  SPFP  to  estimate  the  initial  state  vector 
from  selected  data  points  is  a  gaussian  iterative  scheme  that  requires  two 
geocentric  position  vectors  with  associated  times.  The  initialization  pro¬ 
cedure  implemented  consists  of  the  following  four  steps: 

1.  Selection  of  two  data  points  from  the  refined  data  set  O;  that  is, 

O.  and  O. 
k  J 

2.  Computation  of  two  geocentric  position  vectors  from  selected 
data  points;  that  is,  £^.(t^)  and  r^(L) 

3.  Application  of  the  gaussian  iterative  scheme  to  obtain  cor¬ 
responding  inertial  velocity  vector  jf^(t^) 

4.  Utilization  of  closed  form  expressions  for  the  f  and  g  series  to 
compute  the  state  vector  at  tQ. 

i.  DATA  POINT  SELECTION 

Given  the  refined  data  set  O  =  {O.  :  i  =  1,  .  .  .  ,  n}  and  a  non-zero 
integer  value  for  data  point  selection  index  m,  tte  selection  rule  for  choosing 
k  and  j  is  given  in  Table  IV. 


Table  IV.  Data  Point  Selection  Rule 


m  =  1  , 

k  =  1  and  j  =  n  (first  and  last) 

=  2 

k  =  i  and  j  =  l[n/2]  (first  and  midpoint) 

m  >  3  , 

k  =  I  [n/m  ]  and  j  =  I  [n  •  (m  -  i)/m] 
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2.  POSITION  VECTORS  FROM  SELECTED  POINTS 

Inertial  position  vectors  are  computed  from  selected  data  points  and 

O.  by  applying  two  transformations.  If 
1 


=  {L  j  R.»  A^,  E.} 
S  =  {$,k,h} 


for  i  =  k  and  j 


W8  =  (w8,  W8,  W8) 


then  a  unit  station-vehicle  vector  V1  is  computed  where 


cos  E.sin  A. 
1  l 


=  cos  E.cos  A.  for  i  =  k  and  i 

li  J 


sin  E. 

l 


$  XX 

The  T$  transformation  matrix  is  applied  to  V  to  obtain  vector  U  ;  that  is 


ul  =  t|v1  = 


for  i  =  k  and  j 
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where 


T*  = 


-sin  $ 
0 

cos  $ 


cos  $ 
0 

sin  $ 


This  U1  vector  is  multiplied  by  the  slant  range  R^  to  obtain  vector  Q1 ;  that  is 


Q1  =  R.U1  = 


Q- 


Q. 


for  i  =  k  and  j 


X  X  s  s 

Vector  W  is  computed  by  augmenting  Q  with  (Wj,  0,  W3);  that  is 


wi 

Qi  +  Wi 

II 

°2 

Q3  +  W3_ 

for  i  =  k  and  j 


Finally,  the  inertial  position  vectors  are  obtained  by  applying  the  T® 
transformation  matrix  to  Wl;  that  is 


for  i  =  k  and  j 
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whore 


cos  a 
sin  at 
0 


-sin  at  0 

cos  at  0 

0  i 


and 


at  =  at  +  X  +  t.  u 

So  1  e 


Hence  the  set  of  boundary  conditioxis  {r^,  r.,  At}  is  determined  where 

3 

r^  is  the  geocentric  position  vector  at  observation  time 

r.  is  the  geocentric  position  vector  at  observation  time  t. 

3  J 

til  til 

At  is  the  time  difference  between  the  j  and  k  data  points 
(At  stj  mtk* 

3.  GAUSSIAN  ITERATIVE  SCHEME  TO  OBTAIN  VELOCITY 
VECTOR 

To  determine  a  (two -body)  orbit  satisfying  the  {r^,  r^,  At}  bounding 
'conditions,  the  following  gaussian iterative  scheme  (Ref.  5)  is  used  in  the  SPFP  to 
compute  the  inertial  velocity  vector  r^  from  f  ang  g  closed  form  expressions. 
Compute 


w  =  — J — —  = 

!~j  X  ^kl 


w. 


X 


w. 


w_ 


(cont. ) 


1 


T=n1/2(t.  -y  ^jx1/2At 
r  =  (r  ♦  r  ^2 

j  -j  ~j 


rk  *  (^k  *  ^k} 


cc9<vj  “  vk}  =  (ik  '  — j)/(rkrj)  ’  0  <  (v.  -  vk)  <  7r 

(  x.  • V.  -  x.v.  ) 

sin(v.  -  v,  )  - - - - •  [i  -  cos2(v,  -  V  J]1^2 

J  k  lxkyj  -  -Yki  j  k,J  (22) 

where  (v.  -  v, )  is  the  true  anomaly  difference  between  observation  times  t, 

J  _  K  k 

and  t.,  which  corresponds  to  the  in-plane  angle  difference  for  At. 

If  Wz  5  0,  then  motion  is  direct;  and  if  Wz  <  0,  motion  is  retrograde. 

Obtain  constants  l  and  m  by 

y  4*  r 

^ «  1  A 


m  -  t 


'[2<rj  ’  rk,1/2cos(-J-T^-} 


Note  that  cos[(X  »  Y)/2]  =  +[(1  +  coa(X  -  Y))/2]1/2. 
Choose  y^  =  1,  and  then  compute 


m  „ 
X{i)  =  2  ~  1 

y(i) 


/E.  -  E, 


!  -  2x(i) 
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(E.  -  E, 


)m  ^  [4XU> 


(l-X(i))]1/2 


Equations  (25)  and  (26)  determine  (E.  -  E,  )...  uniquely. 

3  ^ 

Compute  x  from  the  (E.  -  E^)^.  ^  value  from 


(E.  -  E.  )...  -  sin(E.  -  E. ),.. 
■_J . k(l) _  3  k;(i) 

o  /E.  -  E,  \ 


to  obtain  a  new  value  for  y,. , . that  is 

Mx+i) 


y  (i+i )  = 1  +  x(1  +  x(t)' 


The  iterative  cycle  (Eqs.  (24)  through  (28))is  repeated  until 
|Ay|  =  Jy^-+jj  -  ^(i)  I  ~  some  Predetermined  e. 

Computation  continues  with 


1/2  _J,  ,  xi  /2  /Vi  “  Vk\  .  f 

i  1  —  T  2y(rfc  •  r^)  CuSV~J2 - )  smi' 


E.  -  E, 


f  =  i  -f-  [1  -  cos(E.  -  E  )] 

j  3 


g  =  T  -  a3/2[(E.  -  Ek)  -  sin(E.  -  Ek)] 


£k  =  ^1/2<£j  "  fIk)/g 


Hence  ^  and  £k  are  known  and  the  state  vector  is  determined  at  time  tj^. 
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4.  STATE  VECTOR  AT  t 

o 

Closed  form  expressions  for  the  f  and  g  series  are  used  to  evaluate 
the  state  vector  at  I(t^)  =  t^.  If  the  position  and  velocity  at  are  given, 
thfen  the  state  vector  at  t  can  be  expressed  as 


r  =  f  r,  + 
— o  — k 


(33) 


and- 


(34) 


where 


=(f), 


g  r0  sin* 
k 


*  =  v  -  v, 
ok 


f  =  ^ry(*0C0S  $  -  r^sin  $$)  -  (rQsin  $  +  rQcos  $$) 


g  =  (§ )k^0sin  $  +  roCOS  *•) 


$  =  2* 


_h 

2 
r . 


h=  l^kx4 
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b.  TRAJECTORY  GENERATOR 

Trajectory  generation  in  the  SPFP  is  accomplished  by  numerically 
integrating  the  equations  of  motion  and  associated  variational  equations  in  a 
Cowell  formulation  with  an  eighth-order  differencing  me.hod  (Ref.  6).  The 
poeition  and  velocity  components  of  the  state  vector  and  its  associated  partial 
derivatives  are  interpolated  at  each  observation  time  t^  with  a  fifth-degree 
polynomial  -technique . 

i.  EQUATIONS  OF  MOTION 

The  motion  of  the  vehicle  with  respect  to  time  is  described  by  the 
second -order  vector  differential  equation 

£  =  -  “3  +  E  (35) 

r 

with  initial  values 

r(t  )  =  r  ,  r(t  )  =  i 

—  o  — O  —  O  — Q 


and  where 

r  is  a  3-vector  of  rectangular  components  of  position  (x,  y,  s)  in  an 
ECI  reference  system 

r  is  the  magnitude  of  r_  (r  =  |rj) 

[ i  is  the  product  (GM)  of  the  Newtonian  gravitational  constant 

and  mass  of  the  earth 

F_  is  the  acceleration  vector  resulting  from  perturbing  forces; 
i.  e. ,  all  forces  other  than  the  inverse  square  central  force 
due  to  gravitation. 

In  particular,  the  perturbative  acceleration  as  modeled  in  the  SPFP  is  the 
sum  of  the  gravitational  acceleration  due  the  noncentral  force  field  of  the 
earth  F^  and  the  acceleration  due  to  atmospheric  drag  Fg. 
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2.  PARTIAL  DERIVATIVES 

Parameters  associated  v/ith  £  may  be  functionally  '-.elated  to  the  state 
vector  by  the  following  form  of  Eq.  (35) 


I  =— 3-  +Z(£»£»X»t) 

r 


or 


where  y  is  a  vector  of  parameters. 

If  r^  denotes  dr_/dy,  this  leads  to  the  general  vector  variational  equation 

form 


9F 

VW 


with  the  necessary  initial  conditions  r^(tQ),  r^(to),  depending  on  the  parameter 
■y  of  y_.  Numerical  integration  of  this  second-order  vector  variational  equation 
yields  the  partial  derivatives  r^  and  r^  for  each  y  of  £  that  is  a  function  of 
time  and  subject  to  its  respective  differential  equation. 


E.  RESIDUAL  AND  PARTIAL  DERIVATIVE  EVALUATOR 

Once  the  trajectory  point  at  observation  time  t.  is  determined,  then 
data  residual  and  data  partial  derivatives  can  be  evaluated  with  respect  to 
trajectory  parameters.  To  accomplish  this  the  state  vector  at  time  t.  is 
transformed  into  a  data  reference  frame  by  the  transformations  given  in 
Section  III.C.4,  which  allows  an  observation  set  of  data  measurements  to  be 
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computed.  Let  O P  denote  this  set;  that  is 


Op  =  It.,  Ftf,  Ap,  eP, Rp 

i  liiiii 


th 


(36) 


where  the  computed  measurements  for  the  1  observation  are 


Rp  is  the  computed  slant  range 


Ap  is  the  computed  azimuth 


E.  is  the  computed  elevation 

*  Q 

R.  is  the  computed  range  rate 
1.  RESIDUAL  COMPUTATION 


A  residual  set  A.  =  O.  -  O.  is  computed  from  the  refined  observed  data 

Ill  c 

point  Oj  and  the  computed  data  point  O^  at  time  L ;  that  is, 


A.  =\  AR.,AA.,AE.,AR. 

X  1  1  1  1 1 


(37) 


where 

ARj  =  Rt  -  Rp 

AA.  =  A.  -  AP 
x  i  i 

AE.  =  E.  -  Ep 

xii 
•  •  # 

AR.  =  R.  -  Rp 

ill 

' '  v 

These  residuals  of  A^  are  an  index  of  observation  measurement  errors,  force 

model  errors,  and  vehicle  initial  state  vector  errors  from  observation  to 
observation  during  the  data  pass.  Ultimately,  the  differential  correction 
process  will  depend  onA.  for  i  =  1,  .  .  .  ,  n. 
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2. 


DATA  PARTIAL  DERIVATIVES 


Partial  derivatives  of  each  data  measurement  (R,  A,  E,  R)  with  respect 
to  trajectory  parameters  are  denoted  in  the  chain  rule  expression 

90?  90?  9r  90?  9r 

W  “  ~9l  +  ~dT~  h.  ’ 

and  are  necessary  for  the  differential  correction  process  at  each  time  t. 

(i  =  1,  .  .  .  ,  n).  The  sequence  of  computations  leading  to  these  data  partial 
derivatives  is 

a.  Evaluation  of  partial  derivatives  of  computed  observation  set  O1 
with  respect  to  the  state  vector  at  time  ti;  that  is,  90?/9r_  and  c 
90?  /9 r  (refer  to  partial  derivative  equations  in  Section  III.  E). 

£ 

b.  Evaluation  of  90^  /9^_  accomplished  by  using  trajectory  partials, 
9r^/9^  and  9f/9y_,  and  the  observation  partials,  90?/9_r  and 
90P/9r,  to  evaluate  the  chain  rule  expression. 

£ 

The  exact  equations  used  to  evaluate  these  data  partials  90.  /9^  are  given 
in  Section  III.  E. 

F.  RESIDUAL  EDITOR  AND  MATRIX  ACCUMULATOR 

Before  the  accumulation  of  observation  partials  and  residuals  in  matrix 
form,  the  residuals  are  subject  to  the  following  editing  technique. 

1 .  RESIDUAL  EDITING 

As  the  residuals  for  each  observation  are  evaluated,  their  squares 
are  accumulated.  At  the  termination  of  the  accumulation,  the  root  mean  square 
(rms)  is  calculated  for  each  measurement  type  and  used  subsequently  in  the  next 
accumulation  to  remove  (edit)  the  data  points  whose  residuals  have  exceeded 
a  predetermined  multiple  of  the  rms  value.  The  multiple  rms  factor  is 
reduced  for  each  iteration  to  improve  convergence;  hence,  only  the  best 
observations  contribute  to  the  matrix  accumulation. 
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2. 


MATRIX  ACCUMULATION 


Observation  parHals  and  residuals  A.  at  time  t.  are  stored  for 

each  data  measurement  in  a  row  of  the  accumulation  matrix  A«  A  typical  row 
of  A  would  be  of  the  form 


* 


Row  j  (i 


£  j 


*4) 


q  «  th 

where  o.  is  the  computed  j  measurement  for  which 
J 

j  =  1.  range 
j  =  2,  azimuth 
j  =  3,  elevation 
j  =  4,  range  rate 


and 


til 

is  the  h  parameter  with  k  =  i,  .  .  .  ,  7 

til 

A.  is  the  residual  of  the  j  measurement. 


T 

The  weighted  normalized  matrix  product  denoted  by  A  A  is  accumulated  from 

the  A  matrix  at  each  data  point  and  a  diagonal  matrix  W,  which  is  used  to 

normalize  the  units  and  give  relative  weights  to  each  data  measurement. 

T 

The  j,k  element  of  the  A  A  is  expressed  as 


(A  A)jk  =  SAf  jWff  Afk 

f=i 

a°f  i  do£  J  =  1  ’  * r ' 

=h ^ ^  k=i, ... 


7 

7 
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til 

where  o.  is  the  i  data  measurement, 
i 

representatives  of  parameter  vector  x> 
til 

i  data  measurement. 

Note  that  the  elements  of  the  last 
the  form 


Yj  and  y^  are  two  independent 
and  0-j  is  the  weighting  factor  of  the 

T 

column  (or  row)  of  the  A  A  are  of 


(ATA).8 


4 

■E 

i  =i 


S 


T 

which  is  denoted  by  A  B  and  contains  the  nonhomogenous  values  of  the 
linear  system  solved  in  the  differential  correction  process. 

G.  DIFFERENTIAL  CORRECTOR 

Differential  corrections  Ax  to  the  parameters  x  are  obtained  from  the 
solution  of  the  linear  system  described  in  Ref.  1  as 


[ATA  +  ZGTG]Ay.  =  atb 


(38) 


where 

T  T 

A  A  is  the  weighted  normalized  symmetric  matrix  of  the  form  A  WA, 

with  W  denoting  t'he  diagonal  weighting  matrix  and  A  being  the 

matrix  of  observation  partials  and  residuals 

Z  is  an  optimizing  coefficient  to  be  solved  for 

T  T  2 

G  G  is  a  diagonal  matrix  product  of  the  form  (G  G)..  =  g. 

(j  =  1,  ...»  7)  where  G  is  a  diagonal  matrix  ^  ^ 

which  specifies  the  upper  bounds  g.  to  the  corrections 

Ay*(j  =  1,  .  •  .  .  7) 

J 

Ax  is  the  vector  of  parameter  corrections  to  be  solved  for 
B  is  a  vector  of  weighted  residuals 


The  solution  of  this  linear  system,  namely  A^>  is  subject  to  the  constraint 
that  the  optimizing  coefficient  Z  has  a  value  such  that  the  corrections 

=A^(Z)  satisfy  i[GA^_]  [GAyJ  -  l|  <  6,  where  5  is  a  predetermined  con¬ 
stant.  An  iterative  scheme  is  utilized  to  find  suchan  optimizing  coefficient,  but 
is  terminated  after  a  specified  number  of  iterations  if  an  optimum  value  is 
not  found.  The  matrix  of  bounds  G  is  initialized  for  the  iterative  cycle.  The 
modified  G  matrix  depends  upon  the  previously  obtained  value  of  [B  B] ; 
that  is 

if[B^B],<  [B^B],  then  G  =  [m  *A}j] '  (reduction) 


if  [B^B]'  >  [BTB],  then  G  =  [mAy]'  (expansion) 


where 


BXB  is  the  dot-product  of  the  vector  of  weighted  residuals 
(the  sum  of  squares) 

G  is  the  diagonal  matrix  of  bounds 

m  is  a  prespecified  constant 

A>£  is  the  correction  vector 

(prime)^eno^es  va^ue  (or  vecf°r)  obtained  from  the  previous  iteration. 

1.  CONVERGENCE  OF  THE  DIFFERENTIAL  CORRECTION 
PROCESS 

When  a  correction  vector  A^(Z)  is  obtained  in  the  iterative  scheme 
described  above,  then  examination  is  made  of  the  convergence  condition 


ETB  -  IJAA^  -  B  || 2 


where  B  ^  £  and  A^  *  0,  empirically  (that  is,  linear  techniques  are  used  to 
solve  a  problem  that  is  inherently  nonlinear;  therefore,  neither  residuals 
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nor  differential  corrections  are  identically  zero),  and  €  is  a  prespecified 
relative  convergence  criterion..  If  this  inequality  is  satisfied,  the  process 
has  converged;  if  not,  the  iterative  differential  correction  loop  is  continued. 

Terms  associated  with  convergence  levels  are  "current  rms"  and 
"predicted  rms,"  defined  numerically  as 


Current  rms  =  [  (B^Bj/n]  *  ^ 


Predicted  rms  =  (||AA^  -  B||^/n)*^ 


2.  SCHEMATIC  REPRESENTATION  OF  THE 

DIFFERENTIAL  CORRECTION  PROCESS 

The  iterative  differential  correction  loop  depends  on  the  differential 
corrector  that  performs  the  least  squares  process  to  obtain  a  set  of  cor¬ 
rections  A ^  for  the  parameter  y_.  Figure  7  is  a  schematic  representation 
of  the  differential  process  with  emphasis  on  the  function  of  the  differential 
corrector.  The  fundamental  steps  that  identify  the  iterative  differential 
correction  loop  can  be  summarized  as  follows: 

Step  1.  Establish  initial  values  for  the  state  vector  S  (r,  r) 
and  parameter  vector 

Step  2.  Generate  the  vehicle's  trk§^£tory  and  associated  partial 
derivatives  from  S  (r,  r )  and  V. 

O  “ *  *  o 

Step  3.  Evaluate  observation  residuals  A  and  partials  90  /9y 

Step  4.  Perform  residual  editing  and  matrix  accumulation 

Step  5.  Solve  for  the  differential  correction  vector  Aj^  and  test 

for  convergence  "* 

Step  6.  Apply  correction;  that  is,  reset  S  (r,  i)  and  v  -  ^  +  Aj£ 
and  continue  with  Step  2. 
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SECTION  V 


PROGRAM  USAGE 


The  SPFP  has  features  that  permit  the  handling  of  data  in  several  modes. 
Specific  program  control  cards  enable  the  user  to  edit  otherwise  inaccessible 
data,  such  as  that  stored  on  tape.  Program  deck  setup  depends  upon  the  use 
of  the  control  cards.  The  input  or  change  of  constants  is  accomplished  by  a 
generalized  input  routine. 

A.  DECK  STRUCTURE  AND  FORMATS 

The  program  deck  setup  depends  upon  the  types  of  run  desired  and 
requires  consideration  of  station  location  input,  input  of  flags  and  constants, 
observation  data  flocks,  and  program  control  cards.  .  Formats  for  each  of 
these  items  are  described  in  FORTRAN  format  notation. 

1 .  STATION  LO CATIONS 

Up  to  10  stations  can  be  input  to  the  program.  The  station  cards  are 
the  first  items  read  on  any  run  and  contain  identifier  and  station  coordinates 
(latitude,  longitude,  and  height).  The  FORTRAN  format  is  specified  as  (A3, 
10X,  3E15.  8).  A  station  called  TCb  terminates  the  station  card  read  (where 
b  indicates  a  blank  on  the  data  card). 

2.  PRELIMINARY  INPUT 

Following  the  station  location  cards,  the  program  reads  initial  inputs, 
which  include  all  reference  dates,  flags,  and  constants.  The  program  presets 
all  flags  and  constants  to  perform  a  typical  run  with  all  observation- dependent 
constants  in  the  basic  unit  system  of  feet,  degrees,  and  feet  per  second.  The 
epoch  date  corresponding  to  the  earliest  observation  time  is  a  recommended 
input.  The  format  for  this  type  of  input  consists  of  three  fields,  and  as  many 
as  three  sets  of  fields  may  be  placed  on  one  card.  Each  variable  is  identified 
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ff 


by  a  mnemonic  of  up  to  five  characters.  Mnemonics  can  also  designate 
arrays,  subsequent  locations  of  which  can  be  referenced  by  FORTRAN-type 
subscripts.  The  mnemonic  and  the  values  of  the  variable  are  left  justified  on 
the  input  load  sheet.  A  one-character  code  is  used  to  indicate  the  type  of  vari¬ 
able  —  b  for  real  numbers  (floating  point)  and  I  for  integers  (fixed  point).  The 
code  E  terminates  data  input  reading.  The  FORTRAN  format  for  the  value  is 
either  E  or  F  for  real  variables,  and  I  for  (left  justified)  integer  variables. 

A  detailed  list  of  input  mnemonics,  their  preset  values,  and  descriptions  is 
presented  in  Section  V.  1.  B. 

3.  OBSERVATION  DATA  INPUT 

Observation  data  obtained  from  tracking  a  vehicle  can  be  in  either  of  two 
formats,  corresponding  to  those  used  in  the  Aerospace  Generalized  Orbit 
Determination  Program,  TRACE,  versions  D  (Ref.  1)  and  66  (in  preparation). 
Both  formats  contain  station  identification  characters:  the  time  the  observa¬ 
tion  was  measured  (month,  day,  hour,  minute,  and  second);  observation  type 
code;  and  three  data  fields  containing  R,  A,  and  E  measurements  (type  1), 
or  the  first  field  containing  an  R  measurement  (type  7).  The  FORTRAN 
format  for  a  TRACE-D  data  card  is  (A2,  2X,  4F3.  0,  F7.  5,  Ii,  3Fi5.  8),  and 
for  a  TRACE- 66  data  card  (A3,  7X,  4F2.  0,  F8.  5,  IX,  Ii,  2X,  3E15.  8).  Each 
pass  of  data  requires  a  station  identification  TRb  to  terminate  the  reading  of 
observation  data. 

4.  PROGRAM  CONTROL  CARDS 

The  SPFP  control  cards  permit  the  user  to  control  the  program  to  fit  his 
needs.  There  are  four  options  that  tell  the  program;  to  read  a  new  set  of 
flags  or  constants,  to  skip  passes  in  the  observation  data  input  stream,  to 
delete  observations  from  a  data  pass,  and  to  stop  the  program  after  proces¬ 
sing  the  desired  number  of  cases.  It  is  fundamental  to  the  use  of  these  control 
cards  that  one  card  reside  in  memory  and  that  a  new  card  be  read  immediately 
after  the  resident  control  card  is  executed.  The  resident  control  card  is 
continuously  checked  during  the  input  processing  of  the  data.  When  the  con¬ 
ditions  specified  by  the  control  card  are  mat,  the  next  card  in  the  input  stream 
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is  expected  to  be  another  control  card.  The  FORTRAN  format  for  the 
program  control  card  is  (A  10,  315).  Table  V  lists  the*  four  available  control 
cards  in  the  SPFP.  '*“'’“'"4 


Table  V.  Available  Control  Cards  In  the  SPFP 


Name 

Case 

INPUT 

N 

SKIP 

N 

DELETE 

*  N 

STOP 

N 

Counters 


Description 


Add  new  input  to  the  N**1  case 


Skip  cases  N  through  K  (N  <  K) 


Delete  observations  I  through  J 
in  the  N**1  case 


The  DELETE  and  SKIP  cards  are  intended  only  for  editing  continuous  input 
sources;  i.  e. ,  tape. 

Figures  8  and  9  indicate  the  deck  setup  for  runs  from  data  tapes  and 
data  cards,  and  are  followed  by  sample  listings  for  both. 


SAMPLE  DECK  SETUP  RUN  FROM  CARDS 


Monitor  Control  Cards 
Program  Execution 
(End  of  Record  Marker) 


ST 

36, 

,71 

311. 

65 

176.  5 

TS 

RDAT 

1967. 

2 

8. 

3 

24. 

fJT 

E 

ST0P 

5 

1 

IIT 

0 

ISYSIN5 

ST 

8  24  16  43 

0.  0 

1 

3106728. 0 

2.  198 

2.  185 

ST 

8  24  16  44 

o.  o- 

1 

2034792. 8 

5.629 

13.467 

ST 

8  24  16  45 

0.  0 

1 

1163245.' 6 

31. 678 

53.  629 

ST 

8  24  16  46 

0.  0 

1 

1186924. 8 

137.215 

55. 829 

ST 

8  24  16  47 

0.  0 

1 

2074629. 3 

174. 329 

14.  367 

ST 

TR 

8  24  16  48 

0.  0 

1 

3279876. 3 

177. 329 

1.  173 

(End  of  File  Marker) 


SAMPLE  DECK  SETUP  RUN  FROM  TAPE 


Monitor  Control  Cards 
Program  Execution 
(End  of  Record  Marker) 


STA 

TS 

RDAT 

IJT 

E 

SKIP 

DELETE 

SKIP 

INPUT 

RDAT 

E 

ST0P 


1967 


36.71 


1967. 

2 


5 

9 

14 

17 


20 


6 

1 

16 


2 

IIT 


311.  65 

8. 

1 


10 


26. 


(End  of  .File  Marker) 


176.  5 


3  24. 

ISYSIN5 


IMAXIT8 


B.  INPUT  DATA 


All  required  input  is  nominally  set  to  values  in  the  Basic  unit  system  of 
feet,  degrees,  and  feet  per  second.  In  the  following  list,  each  variable  is 
described  and  the  preset  value  designated  in  the  value  field. 


PRESET  CONSTANTS  AND  FLAGS 


Code  Location 

Value 

Description 

RDAT 

1964. 

Epoch  year 

2 

1. 

Month 

3 

0. 

Day 

REFE 

. 000312 

Elevation  refraction  correction  coefficient 

REFR 

900. 

Range  refraction  correction  coefficient 

SCRNG 

1. 

Range  scaling  factor 

RMIN 

0. 

Minimum  range 

RMAX 

1.E30 

Maximum  range 

EMIN 

0. 

Minimum  elevation 

TPASS 

1440. 

Maximuin  duration  of  pass  (min) 

I 

SYSIN 

5 

Card  reader  input  file 

I 

SYS0U 

6 

Printer  output  file 

I 

JT 

2  (5) 

Observation  data  input  tape  unit  (card  input) 

I 

IT 

1  (0) 

Data  format  TRACE- 66  (-D) 

I 

NDIC 

5 

Initial  condition  determination  selector 

EPS 

.  0001 

Two-body  orbit  iteration  convergence  criterion 

SIGMA 

100. 

Range 

2 

.  1 

Azimuth 

3 

.  1 

Elevation 

4 

0. 

Range  rate  data  weighting 
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Code  Location  Value 

» 


Description 


BNDS 

50000. 

2 

50000. 

3 

50000. 

4 

500. 

5 

500. 

6 

500. 

7 

0. 

SCBNi 

.  5 

SCBN2 

2. 

EDITN 

6. 

SCEDN 

.8 

EPSLN 

.  001 

TR 

50000. 

2 

5. 

3 

5. 

4 

0. 

MAXIT 

5 

JPN 

10 

IPN 

1  (2) 

IW 

1 

DF 

VF 

F 

GM 

XNG 

XMG 

CNM 

SNM 


Correction  constraints 


Dragj 

Reduce  f  Bounds  each  iteration 
Expand 

n  -  <r  rms  residual  editing  factor 
Scale  editing  factor  each  iteration 
Relative  convergence  criterion 
Range  'j 
Azimuth  j 

f  Initial  residual  acceptance  limits 


5.  Azimuth 

-  •  Initial  residual  acceptance  limit! 

5.  Elevation  r 

0.  Range  rate 

5  Maximum  number  of  iterations 

10  Number  of  observations  to  be  generated  (£  10) 

1  (2)  Print  (and  punch)  generated  observations 

1  0  print  iteration  summary  only 

1  print  residuals 

2  print  the  matrix  of  partial  derivatives  and 
residuals  at  each  observation  point 

20925738.  Distance  conversion  factor 

348762.  3  Velocity  conversion  factor 

.  335233E-2  Earth  oblateness  coefficient  (=  i /298.  3) 

.  55303935E-2  Gravitation  constant 

2.  n  ] 

0  n  I 

>  28-term  geopotential  input 
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Location 

Value 

Description 

NT 

r 

Number  of  geopotential  terms 

JN0RM 

i 

Normalization  of  geopotential  coefficients 

HO 

1. 

Initial  integration  step  size 

HMIN 

0.  015625 

Minimum  step  size  (1/64) 

HMAX 

64. 

Maximum  step  size 

ER 

i.E-12 

Integration  truncation  control 

-57- 


REFERENCES 


1.  C.  C.  Tonies,  Trace  Orbit  Determination  Program  Version  D, 
Report  No.  TR-669  (9990)  -3,  Aerospace  Corp.,  El  Segundo, 
California  (September  1966). 

2.  The  U.  S.  Standard  Atmosphere.  National  Aeronautics  and 
Space  Administration,  Washington  D.  C. ,  (1962), 

3.  Samuel  Pines,  "Variation  of  Parameters  for  Elliptic  and  Near 
Circular  Orbits,"  Astronomical  J. ,  66,  No.  1,  (No.  1286), 
(February  1961). 

4.  R.  A.  Minzner,  K.  S.  W.  Champion,  and  H.  L.  Pond,  The  ARDC 
Model  Atmosphere,  1959,  Air  Force  Surveys  in  Geophysics  No.  115 
Air  Force  Cambridge  Res.  Center  (August  1959). 

5.  P.  R.  Escobal,  Methods  of  Orbit  Determination,  John  Wiley  &  Sons, 
Inc.  ,  New  York,  (1965). 

'6  J.  F.  Holt,  ASC  DEQ4  Floating  Point  Cowell  (Second  -  Sum),  Runge 
Kutta  Integration  of  Second  Order  Ordinary  Differential  Equations, 
ATM  67  (9990) -89,  Aerospace  Corp. ,  El  Segundo,  California 
(June  1967). 


$ 

Not  available  outside  Aerospace  Corp. 


-59- 


UNCLASSIFIED _ 

Security  Classification 


DOCUMENT  CONTROL  DATA  -  R&D 


(Security  classification  ot  title  body  of  abstract  end  indexing  annotation  mutt  be  entered  when  the  overall  report  it  classified) 


1  ORIGIN  — TIN  0  ACTIVITY  (Corpotato  author) 

2a  report  security  classification 

Aerospace  Corporation 

Unclassified 

El  Segundo,  California 

Zb  CROUP 

3,  REPORT  TITLE 


SINGLE  PASS  FIT  PROGRAM  (SPFP) 


4  DESCRIPTIVE  NOTES  (Type  ot  report  and  Inclusive  datee) 


5  AUTHORS  (Last  name,  first  name,  initial) 

Buechler,  Guenther  and  Walker,  Dennis  G, 


6-  REPORT  DATE 

7a  total  no.  of  PAGES 

65 

7b.  NO.  OF  REFS 

6 

6a  CONTRACT  OR  GRANT  NO 

F04695-67-C-0158 

6.  PROJECT  NO. 

c 

d 

9a.  ORIGINATOR'S  REPORT  NUMSERfSj 

TR -0158(31 10-01  )-l 

$b.  OTHER  REPORT  NOfS>  (A  ny  other  numbers  that  may  be  as  tinned 
this  report) 

SAMSO-TR-68-22 

10  AVAIL  ARILITY/LIMITATIQN  NOTICES 

This  document  is  subject  to  special  export  controls  and  each  transmittal  to 
foreign  governments  or  foreign  nationals  may  be  made  only  with  prior  approval 
of  SAMSO(SMON) _ 


11  SUPPLEMENTARY  NOTES 

12.  SPONSORING  MILITARY  ACTIVITY 

Space  Systems  Division 

Air  Force  Systems  Command 

Los  Angeles,  California 

1  13-  ^ABSTRACT  1 

The  Single  Pass  Fit  Program  (SPFP)  is  designed  for  use  on  the  CDC  6400  or 
6600  machine  as  the  principal  computational  tool  for  preliminary  orbit  deter¬ 
mination  and  processing  of  single  arcs  of  raw  tracking  data.  By  an  iterative 
differential  correction  procedure,  it  solves  for  a  set  of  trajectory  parameters, 
namely,  initial  position  and  velocity  (and  drag  coefficient  if  desired)  of  the 
vehicle,  which  minimize  the  differences  between  measured  and  computed  obser¬ 
vation  measurements.  A  comprehensive  description  of  the  mathematical  model 
implemented  in  the  .SPFP  is  presented  with  emphasis  placed  on  key  character¬ 
istics,  such  as  reference  systems  and  equations,  data  processing  and  self- 
initializing  techniques,  residual  editing,  and  differential  correction.  In  addi¬ 
tion,  a  basic  usage  guide  is  given  complete  instructions  for  input  data  prepar¬ 
ation  and  program  operation. 


DO  '-0RM 

(FACSIMILE! 


1473 


UNCLASSIFIED 


Security  Classification 


KEY  words 


rr«: 


Differential  correction 
Least  squares 
Orbital  Determination 
TRACE 

Trajectory  generation 
Single  pass  fit 

Preliminary  orbit  determination 


Abstract  (Continued) 


